o 

(N 

> 
O 



43 

I 

-I— > 

cr 



> 

m 
o 

(N 

On 
O 
O 



Efficient Algorithm for Asymptotics- Based Configuration-Interaction Methods and Electronic Structure of 
Transition Metal Atoms 
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Asymptotics-based configuration-interaction (CI) methods [G. Friesecke and B. D. Goddard, Multiscale Model. 
Simul. 7, 1876 (2009)] are a class of CI methods for atoms which reproduce, at fixed finite subspace dimension, 
the exact Schrodinger eigenstates in the limit of fixed electron number and large nuclear charge. Here we develop, 
implement, and apply to 3d transition metal atoms an efficient and accurate algorithm for asymptotics-based CI. 

Efficiency gains come from exact (symbolic) decomposition of the CI space into irreducible symmetry subspaces at 
essentially linear computational cost in the number of radial subshells with fixed angular momentum, use of reduced 
density matrices in order to avoid having to store wavefunctions, and use of Slater-type orbitals (STO's). The required 
Coulomb integrals for STO's are evaluated in closed form, with the help of Hankel matrices, Fourier analysis, and 
residue calculus. 

Applications to 3d transition metal atoms are in good agreement with experimental data. In particular we reproduce 
the anomalous magnetic moment and orbital filling of Chromium in the otherwise regular series Ca, Sc, Ti, V, Cr. 

PACS numbers: 31.15.ve, 31.15.vj, 02.70.Wz, 31.15.-p, 32.30.-r 



I. INTRODUCTION 

The search for accurate computational methods for 
the iV-electron Schrodinger equation at moderate com- 
putational cost has been a focus of activity for several 
decades- -5 . The present article is a contribution to one 
part of the picture, wavefunction methods for atoms. We 
develop, implement, and apply to transition metal atoms 
an algorithmic framework which renders asymptotics- 
based Configuration-Interaction (CI) computations for 
atoms with basis sets of up to 50 one-electron spin or- 
bitals, up to 30 electrons, and full resolution of all va- 
lence electron correlations feasible. An attractive feature 
of our framework is that many steps are done symboli- 
cally, by building upon, systematizing, and automatizing 
the paper-and-pencil analysis of asymptotics based CI 
for small atoms and minimal bases in Ref. [J A Mat- 
lab/Mathematica implementation is available at Ref. 0. 

CI methods^ approximate the electronic Schrodinger 
equation by projecting it onto a well chosen subspace 
spanned by Slater determinants. More precisely, the 
Schrodinger equation for an atom or ion with N elec- 
trons is 

H^ = Et/j, ^ € L^(M 3 x{-i ) +i}) JV ) (1) 

where H is the Hamiltonian of the system, see © below, 
ij) the wavefunction and E the energy. The wavefunction 
ip = ip(xi,si, . . . ,xn,sn) depends on the positions Xi € 
R 3 and spins s, G {— |, +|} of all electrons, and belongs 
to the space L^((R 3 x {— |, +^}) N ) of square-integrable, 
antisymmetric functions on (K 3 x {— 5>+§}) • A CI 
method is an approximation of (fj]) by an equation of 
form 



PHPV = Et>, t/j € V C L 2 a ((m 3 x {-!,+!}) 
P = orthogonal projector onto V, 



N 



(2) 



where V is a Span of a finite number of Slater determi- 
nants \xij ■ ■ ■ Xi N ) built from a finite number of spin or- 
bitals {xi, • • ■ , Xk} C L 2 (R 3 x {-±, +i}). We recall the 
well known fundamental difficulty of CI methods: Eq. ([T]) 
is a partial differential equation in very high space di- 
mension, e.g. dimension 72 in case of a single Chromium 
atom as treated in this paper. Hence when discretizing 
the single-electron state space by a reasonable number 
of spin orbitals, L 2 (R 3 x {— |, +|}) « Span{xi, -,Xk}-, 
the ensuing natural choice V — Span {\xh 1 1 1 Xin) '■ 1 — 
i\ < ■ ■ ■ < ijv < K} (full CI) has a prohibitively large 
dimension, (^J . 

Our principal contribution here is the development of 
an efficient algorithm that minimizes the curse of di- 
mension. The main savings come from exact (i.e. sym- 
bolic) and efficiently automated exploitation of symme- 
try to perform dimension reduction. Other ingredients 
are use of reduced density matrices in order to avoid 
having to store wavefunctions, and use of Slater- type or- 
bitals (STO's) including exact orthonormalization and 
Coulomb integral evaluation. The algorithm has been 
implemented for a recent variant of CI, asymptotics- 
based Cli., which exploits the asymptotic results in Ref. 8 
and has the attractive features that the CI subspace, 
if its dimension is K, reproduces correctly the first K 
Schrodinger eigenstates in the limit of fixed K, fixed elec- 
tron number, and large nuclear charge Z. (This limit, 
which has a large literature (see in particular Ref. M 
and ITol ) captures the physical environment of inner shell 
electrons, and has the multiscale property that the ra- 
tio of first spectral gap to ground state energy of the 
Schrodinger equation tends to zero^, with the experimen- 
tal ratio for true atoms being very close to zero, about 
1 part in 1000 for Carbon and Oxygen and 1 part in 30 
000 for Cr and Fe.) The main part of the algorithm, au- 
tomated symmetry reduction, can be easily adapted to 
other CI methods and orbitals (such as Gaussians). 

As a typical application, we treat here the transition 
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metal series Ca, Sc, Ti, V, Cr, modelled by 18 core elec- 
trons occupying Slater orbitals of type Is to 3p, and an 
active space consisting of 3d, 4s, 4p, Ad Slater orbitals 
(of either spin) accommodating the two to six valence 
electrons. The resulting CI space V for Cr has dimen- 
sion d — ( 2 6 8 ) = 376740, and the CI Hamiltonian has 
d(d + l)/2 7 x 10 10 entries. But automated symme- 
try reduction shows (see Table IIIII) that only 14 basis 
functions contribute to the experimental ground state 
configuration and symmetry, [Ar]4s 1 3c? 5 7 S, allowing to 
evaluate the ensuing eigenvalues and -states easily and to 
machine precision. Our results, detailed in Section IVIII 
below, provide an ab initio explanation of the anoma- 
lous magnetic moment of Chromium (experimentally, the 
ground state has six instead of the expected four aligned 
spins) and the underlying anomaly in the filling order 
of 3d versus 4s orbitals in the semi-empirical orbital pic- 
ture of transition metal atoms (Chromium, unlike its four 
predecessors Ca, Sc, Ti, V, possesses only one instead 
of two 4s electrons). It is well knownii~— that single- 
determinant Hartree-Fock, relativistic Hartree-Fock, and 
density functional theory calculations (even with the best 
exchange-correlation functionals such as B3LYP) render 
the correct filling orders and ground state symmetries 
only for some but not all transition metal elements (see 
Section EED- 

In the remainder of the Introduction we describe 
our algorithm for exact (symbolic), efficient symmetry 
partitioning. The (non-relativistic, Born-Oppenheimer) 
Hamiltonian 

i=l v ' ' l<i<j<N 1 Jl 

governing atoms/ions with N electrons and nuclear 
charge Z has the symmetry group 

SU(2) x 50(3) x Z 2) (4) 

consisting of simultaneous rotation of electron spins, and 
simultaneous rotation and sign reversal of electron posi- 
tions. This leads to the well known conservation law that 
the Hamiltonian leaves the simultaneous eigenspaces of 
the spin, angular momentum and parity operators 

L 2 , L„ S 2 , S z , R (5) 

invariant (see Section III Al for precise definitions of these 
operators). The fact that partitioning into symmetry 
subspaces significantly lowers computational costs has 
long been known to, and exploited by, theorists (see e.g. 
Ref. [TEh . A striking example is the paper-and-pencil sym- 
metry decomposition^^ of a minimal asymptotics-based 
CI Hamiltonian PHP for the second period atoms He to 
Ne, with active space consisting of the eight 2s and 2p 
spin-orbitals accommodating the valence electrons. For 
Carbon, there are four valence electrons, so the active 
space has dimension (,) = 70, and the CI Hamiltonian 



is a 70 x 70 matrix. But due to symmetry it decomposes 
into fifteen 2x2 blocks and fourty lxl blocks. 

The main algorithmic steps which automate such de- 
compositions are as follows. 

(a) One starts by partitioning the CI space into config- 
urations, i.e., subspaces like ls™ 1 2s™ 2 2j/ 13 . . . with a fixed 
number of electrons in each subshell (see Section III Bl 
below). It suffices to symmetry-decompose each config- 
uration, because the symmetry group, unlike the Hamil- 
tonian, leaves each configuration invariant individually. 

(b) Each configuration is isomorphic to a non- 
antisymmetrized tensor product of lower-dimensional fac- 
tors. The tensor factors consist of single Is, 2s, 2p, . . . 
subshclls. See Section III Bl This product structure is 
essential for Step (d) below. 

(c) The splitting up of each factor into simultane- 
ous eigenspaces of the symmetry operators ([5]) is done 
via a suitable algorithm from the mathematics litera- 
ture for simultaneous diagonalization of commuting ma- 
trices, for instance that of Bunse-Gerstnert, Byers and 
Mehrmann 17 . (We are indebted to Folkmar Bornemann 
for helpful advice regarding this step.) Exact eigenstates 
are recovered from the numerical eigenstates through ex- 
ploiting that the squares of the eigenstate coefficients are, 
by representation theory, rational numbers. 

(d) Given simultaneous eigenstates of the symmetry 
operators for each factor, simultaneous eigenstates of a 
two-factor tensor product are known explicitly in terms 
of the well-known Clebsch-Gordan coefficients, and those 
for a many-factor tensor product are easily obtained by 
iteration of the Clebsch-Gordan formulae. This yields 
the desired decomposition of each configuration. 

A key feature of the algorithm (a) , (b) , (c) , (d) is the 
computational cost grows only linearly with the number 
of subshells, provided the angular momentum cutoff is 
held fixed. Thus, say, the cost of including s orbitals of 
type Is, 2s, . . . , ns is only 0(n). See Section IVTl 

The structure of this paper is as follows. In Section [TT1 
we briefly review asymptotics-based CI. Section |HT] con- 
tains the main contribution of this paper, namely exact 
reduction steps leading to significant savings of compu- 
tational time and memory storage. In Section IIVI treat 
orthonormalization and Coulomb integral evaluation for 
general atomic Slater-type orbitals. We summarize all al- 
gorithmic steps in Section |Vj and carefully estimate the 
costs in Section IVII Finally, in the last section we apply 
the algorithmic framework to the electronic structure of 
potassium, calcium and the transition metals scandium 
to zinc. 



II. ASYMPTOTICS-BASED CI 

We briefly recall the set-up and features relevant to 
the present work, referring to Refs. [l] and[l6| for further 
information. 
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A. Symmetries 

Due to invariance of the Hamiltonian ([3]) under the 
symmetry group (j4} , the set of operators (JS]) commutes 
with the Hamiltonian and with each other, for arbitrary 
N and Z . These operators play an important role in our 
algorithmic framework. Here and below we use the stan- 
dard notation L = X)t=i ^(*) (many-body angular mo- 
mentum operator), L(i) — Xi A AV Xi (angular momen- 
tum operator acting on the position coordinates Xi G M 3 
of the iih electron), L x , L y , L z (components of L), and 
analogously for spin (see e.g. Ref. [la). The parity op- 
erator ip(x\, st,..., a; at, sn) h> ip(—xt, St, ■ ■ ■ , —xn, sjv) 
is denoted by R. 



B. Configurations 

Our treatment of symmetry reduction is independent 
of the particular orbitals used, and works within the con- 
text of general iV-electron configurations as introduced in 
Ref. [H, definition 2.2: Let 

V u V 2 ,---cL 2 (R 3 x {-£,+§}) (6) 

be any collection of mutually orthogonal subspaces of 
the single-electron Hilbert space, which are irreducible 
representation spaces for the joint spin and angular mo- 
mentum algebra Span {L x , L y • L z , S x , S yi S z\ (or, equiv- 
alently, which are joint eigenspaces of L 2 and S 2 with 
minimal dimension (2£ + l)(2s + 1) given the respective 
eigenvalues £(£+1) and s(s+l)). Then, a configuration of 
an iV-electron atom or ion is a subspace of the antisym- 
metrized A-electron state space L 2 ((^ 3 x {-^,+^}) N ) 
of the following form: 

C d u ...,d k = Span{|xi,...,Xiv) : {Xi,---,Xn} any ON 
set in L 2 (K 3 x {-\,+\}) with j) {i : X i e Vj} = dj} 

where (dt, d/.) — d is a partition of N (i.e. < dj < 
dim Vj , J2j dj = N) . Physically, the Vj are subshells and 
the dj are occupation numbers. The main point is that 
all choices of the Xi' s consistent with the requirement 
that a fixed number of them have to be picked from each 
Vj have to be included. 

Configurations, unlike general subspaces of L 2 ((R 3 x 
{— i +^}) N ) spanned by a basis of Slater determinants, 
are invariant under the symmetry group (jj) and its gen- 
erators L, S and R, and in particular under the oper- 
ators ©. The same holds for multi-configuration sub- 
spaces 

y = Span{c d(1, ,...,C d<M> } (7) 
where each C d< ' is a configuration. 



C. Asymptotics-based selection of 
configurations 

Eq. ([7]) still leaves a great deal of freedom for the pre- 
cise specification of the CI subspace V. In asymptotics- 
based CI 1 , the traditional step of an intermediate 
Hartree-Fock calculation to determine orbitals is replaced 
by the theoretical requirement that the ansatz space re- 
produce correctly the lowest Schrodinger eigenstates in 
the iso-electronic limit Z — > oo (see Theorem Q] below) . 
This 

- requires Slater type orbitals (STO's) instead of the 
asymptotically inexact linear combinations of Gaussians 
which are common in molecular calculations, 

- and corresponds to full CI in an active space for the 
valence electrons (instead of truncating valence electron 
correlations in terms of order of excitation with respect to 
a reference determinant as in double-excited CI, or non- 
linearly approximating them as in coupled cluster the- 
ory). 

Asymptotics-based CI preserves the spin and angular 
momentum symmetries of the original Hamiltonian, and 
obeys the virial theorem, by determining orbital dilation 
parameters self-consistently for the actual CI wavefunc- 
tions instead of precomputing them via a Hartree-Fock 
calculation. (The fact that methods with self-consistent 
dilation parameters always obey the virial theorem was 
pointed out by Lowdini^.) 

The specific asymptotics-based CI model for atoms 
used in this paper is as follows. 

(A) (Choice of a parametrized, asymptotically exact 
family of subspaces) We specify the orbital spaces 
in Eq. ([6]) as 

V^ e := Spaxi{ip n im1;ipTUmi} m= -t...e , 
n = l,2,..., £ = 0,...,n-l 

with orthonormal Slater (or hydrogen-like) orbitals 

i>n£m(x) = r l Y lrn {{}, If) p n i {Z u , Z nt , r) e 

r = |x|, and polynomials p n t(Z u , . . . , Z ni , ■) of 
order n — I — 1, see equation (f2Tj) . Here, Z = 
(Ztfi, Z 2 fi, Z 2i i, . . . ) is a vector of dilation parame- 
ters Z n £ > 0. We then set 

V z := Span I [j C d \ 
ldex> J 

where d = (d n i), n = 1, 2, . . . , I = 0, . . . , n — 1 is 
a vector of occupation numbers which sum to N, 
and I? is a finite set of such vectors such that 

(i) < d^ < dim V% = 2 • (21 + 1). 

Prototypical is the set V consisting of all configu- 
rations such that, with respect to alphabetical or- 
dering of the indices (n,£), 
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(ii) d nt = for (n, £) > (n, £) max 

(iii) dra = 2 • (2t+ 1) for (n,*) < (n,£) min . 

Here (ii) is a cutoff condition, and (iii) says that all 
subshells up to (n, £) m i n are completely filled. 

(B) (Subspace eigenvalue problem) For each symmetry 
subspace 



V z - 

v lsp - 



{4>E V z : L 2 4> = 1(1 + l)i/>, } 



S 2 il) = s(s + 1)"0, i? V = .P "0 } 

of y z (angular momentum, spin and parity quan- 
tum numbers I, s and p, respectively), E\ 
eigenvalues of P z H P z on V? ' , ipf 1 ' 
sponding orthonormal eigenstates, where P z = or- 
thogonal projector of Ll((R 3 x {-|, +^}) N ) onto 
V z . 



ci.z 



corre- 



(C) (Variational parameter determination) For each 
symmetry subspace V z „ C V z , 



argmin z (minj .E- 

, CI,Z, 



CI,Z\ 



v tsp 

Ef 



Here argmin 2 ,/(a;) denotes a minimizer of the functional 
/. We remark that minimizing dilation parameters Z are 
expected to exist provided the nuclear charge Z is greater 
or equal to the number TV of electrons (in which case 
the full Rayleigh-Ritz variational principle possesses a 
minimizer—). In our numerical computations we always 
found this to be the case. 

Also, for future reference we define 

c CI := number of core spin-orbitals of the CI model 
= £ 2(2*+ 1), 



(n,£)<(n,£)„ 



t := total number of spin-orbitals of the CI model 
= Y, 2-(2£+l). 

(",')<(n,f)m« 

Of course, the model only makes sense (i.e., the space V 
is nonempty) provided the cutoffs (n,£) m - m , (n,£) max are 
chosen so that c CI < N < t GI . 

We summarize the asymptotic properties of the above 
model in the following straightforward generalization of 
Theorem 2.1 in Ref. [0 on second-period atoms. The 
following numbers associated with the non-interacting 
TV-electron atom play a role: n_(iV), n+(N), c(N), 
t(N) which denote the number of closed shells, closed or 
open shells, core spin-orbitals, and core or valence spin- 
orbitals, respectively. Explicitly 16 , n_ and n + can be 
expressed in terms of the number of spin-orbitals in the 

first n> hydrogen shells, /(»') := J2to 2 ' ( 2£ + 1). 

as the largest integer such that /(n_) < N, respectively 
the smallest integer such that f(n+) > N. One then has 
c(N)=f(n-(N)),t(N) = f(n+(N)). 



Theorem 1. (Correct asymptotic behaviour) The CI 
model (A), (B), (C) with T> given by (ii), (iii) has the 
following properties. In the large nuclear charge limit 
Z — > oo for N and dim V z fixed, the lowest 



min (t(N), t CI ) - max (c(N), c CI ] 
N -max(c(/V),c CI ) 



eigenvalues E^ 1 < E§ 1 < . . . and E% < E% < . . . (re- 
peated according to multiplicity) of the CI model respec- 
tively the Schrodinger equation ([1]) satisfy 

Ef 1 
lim-^- = 1. 

e j 

If moreover c CI < c(N) (i.e. the CI model does not 
constrain the occupation numbers of any non-core or- 
bitals) and t cl > t(N) (i.e. at least all core and valence 
orbitals are included in the CI model), then there exist 
orthonormal CI respectively Schrodinger eigenstates ip^ 1 
and ipi corresponding to the above eigenvalues such that 

i.Ci 



lim||Vi -M = 0, 
is the norm on the N -electron space L 2 a ( 

N\ 



where 

{-mm 

Finally, under the same condition on c 0i an 
spectral gaps satisfy 

AEf 1 
lim = 1 

AE* 



ci „„ d t ci tht 



(8) 



whenever AEj > 0, where AEf 1 = E^-E^ 1 and AE = 
Ej - E x (j >2). 

We emphasize that the above theorem only covers the 
regime of large Z. For neutral atoms, the highest eigen- 
states in the ansatz space of asymptotics-based CI are 
typically observed to lie above higher Rydberg states or 
even above the bottom of the continuous spectrum. 



III. EXACT REDUCTION STEPS AND LS 
DIAGONALIZATION 

This section explains exact reduction steps which are 
essential for cutting down the calculation time and stor- 
age requirement of the algorithmic implementation. 



A. Tensor product structure of configurations 

Our first observation connects TV-particle configura- 
tions (see Section III Bp to the non-antisymmetrized ten- 
sor product of antisymmetrizcd dj-particle states, pre- 
serving the action of the angular momentum and spin op- 
erators. Here and below, we use the standard notation 20 
A n V for the n-fold antisymmetrized tensor product of a 
vector space V, and V <S> W for the tensor product of two 
spaces V and W. 
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Proposition 2. Consider irreducible representation 
spaces Vl, . • . ,Vk as in equation ([6|) and particle num- 
bers d±,...,dk > 0. Then the following isometric iso- 
morphism holds, 



Qd 1 ,...,d k 



iA dj Vi. 



(9) 



3=1 



A canonical mapping of basis vectors is given by 



T:\xl 



1 Xd\ 5 



\xl, 



\xl 



with xl G Vj f or M hj- Moreover, T commutes with 
the action of the angular momentum and spin operators, 



i.e., 



/or aZi V € C * where L, S on the left hand sides are 
N -particle operators and each Lj , on i/ie rig/ii hand 
side acts on the dj -particle tensor factor A dj Vj ■ 



Proof. Clear from the definitions. 



□ 



In particular, dim(C dl — ><** ) = C""^) • N «te that 
equation @ inherently takes into account the antisym- 
metrization of fermionic wavefunctions, without requir- 
ing any additional normalization factors. The isome- 
try ^ is reflected in the algorithmic implementation 
by ordering Slater determinants lexicographically and ar- 
ranging coefficients accordingly, see Figure [T] 



fci 



K\ + li2 



) 



di 



FIG. 1. Lexicographical ordering of the 



Slater de- 



terminants restricted to a fixed configuration involving two 
one-particle subspaces Vi, Vi of dimension fci, k-z, with d\ 
particles in orbitals 1, . . . ,ki and di particles in the remain- 
ing orbitals fei + 1, . . . , ki + fe. Filled circles correspond to 
filled orbitals. Usefully, restriction to a configuration pre- 
serves lexicographical ordering. 



B. LS diagonalization 

A priori, the diagonalization of the angular momentum 
and spin (LS) operators seems as expensive as diagonal- 
izing the Hamiltonian itself, yet it turns out to come at 
much cheaper costs. It involves mostly algebra, and can 
be done prior to setting up the Hamiltonian. 

1. Calculate all irreducible LS-eigenspaces for each 
many-particle subshell. In more detail, let u — s,p,d, . . . 
denote the angular momentum subshells in common 
chemist's notation and set 

V u := Sp(iTi{Y em t,Y em l} e=angiuhm= _ L i , 

with the spherical harmonics Yi m . Note that this is an 
explicit realization of the spaces in Eq. © . Then, for all 
n= 1, ... ,dim(V u ) (equal to 2x (2£+lj), decompose the 
n-particle space A n V u into the direct sum of irreducible 
spin and angular momentum representation spaces. That 



A n K = 0K- 



(10) 



such that 



L 2 ip 
S 2 (f 



dim(K 



U ik + 1) <p, 

s l (si + 1) if V<p e V u - 

(21, + 1) (2s, + 1) . 



Explicit results are shown in Table |TJ subshells from 
A 4 Vrf to A 10 Vd are omitted for brevity's sake, and only 
states with maximal L z and S z quantum numbers are dis- 
played; applying the ladder operators L_ = L x —iL y and 
S— = S x — iS v yields the remaining wavefunctions. Note 
that symmetry levels can appear twice within a many- 
particle subshell, e.g., 2 D in A 3 Vd- In concordance with 
the Clebsch-Gordan method below, the ordered single- 
particle orbitals are L z and S z eigenstates, denoted by 

(s,s) for V s , 

(pi,Pi,Po, Po, Pm , prf) for V p , 

(d 2 ,d2,di, di, . . . ,dra, dra) for V d . 

The highest quantum number appears first, and ~ equals 
spin down 1 (convention as in Ref.El). 

The decomposition (|I0[) first requires a matrix repre- 
sentation of the angular momentum and spin operators 
L x , L y , L Z ,S X , S y , S z on A n V u . Obtain it by starting from 
the canonical single-particle representation on V u (spher- 
ical harmonics) and writing the ro-body operator in the 
form B = J2i jbij a t a j, where 6y are the coefficients of 
the single-particle representation and a\ , aj are fermionic 
creation and annihilation operators. The operators af a j 
map Slater determinants to Slater determinants; thus all 
entries of their corresponding matrix representation are 
or ±1. 
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The next task to arrive at (|T0|) involves the simultane- 
ous diagonalization of the pairwise commuting operators 
L 2 ,S 2 \L Z ,S Z . We present two alternatives. 

Alternative 1 

apply an algorithm of choice, e.g. Ref. [I?], for the si- 
multaneous diagonalization of L 2 ,S 2 ,L Z , S z on A n V u , 
denoting the eigenvalues or "quantum numbers" of 
L 2 , S 2 ,L Z , S z by £(£ + 1), s(s + 1), me, m s , respectively 

for each subspaces W with me = £ and m s = s do 
choose an ONB {ipi, . . . ,ip r } of W 
for j = 1, . . . , r do 

add V un j := Span{</jj, L-ipj, S-tpj, L^S-ipj, . . . } to 
the decomposition (fT0|) 
end for 
end for 

Note that the iterative application of the ladder operators 
L_ and S- ensures that the resulting subspaces V un j are 
invariant irreducible representation spaces. 

Alternative 2 

count <— 

f I 3 n odd 
for £ = 0,1,... ands= \ do 

[ 0, 1, . . . n even 

X := (L 2 - £{£ + 1) id | S 2 - s(s + 1) id | 
L z - ^id | S^-s id) 

calculate W := Kern (ll*) C A n K 
if W = then 
next for loop 
end if 

choose an ONB {<^i, . . . , tp r } of FT 
for j — 1, . . . ,r do 

add Kin; := Spanjt^j, L^ipj, S-<pj, L-S~tpj, . . . } to 
the decomposition (fTO]) 
end for 

count count + r 

stop if ^ -°" nt dim (Kni) = dim (A n V u ) 
end for 

The second alternative exchanges the direct diagonaliza- 
tion in alternative 1 for testing all potential eigenvalues 
(that is, integer or half-integer numbers) with me = £ 
and m s = s. Efficient numerical methods exist for com- 
puting the kernel W, which take advantage of the sparse 
structure of the matrix representation. 

2. Consider iV-electron configurations C ni '"' ,nk as- 
sembled from the above single-particle subshells, with 
nj electrons in subshell j (angular momentum Uj) such 
that N = X^ n j- Using the decomposition in step 1, 
simultaneously diagonalize the pairwise commuting op- 
erators (O acting on C n i>— > n <=. (The parity operator R is 
constant on C™ 1 '""'™* anyway and needs no further con- 
sideration.) In more detail, the isometry © and the 
decomposition p0[) imply 

£»!,...,»*„ i,. Vl , (g)i; (ii) 

I=(il,...,i k ) j 



Config 


Sym 


L z Sz 




a'K 


9 r-i 

2 S 


„ -i 
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aW p 
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ipi) 
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\piPo) 
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2 po 
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2 i 




A 4 V P 







75 (- PiPiPmPni) - PiPoPoPm) 
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TABLE I. Irreducible LS eigenspace decompositions of A n V u 
in Eq. (|10[1 . showing states with maximal L z and S z quantum 
numbers only. 



By construction, each Vi is uniquely characterized by its 
eigenvalues with respect to the LS-operators L 2 and S 2 
acting on the jth tensor factor. Since 

3 3 



G 



all operators 

L ) S , L Z) S z , R, Lj, Sj j = 1, . . . , k 

commute pairwise, and it follows that each Vj is an in- 
variant subspace of the operators §5$). Thus, the diago- 
nalization can be performed on each Vj independently. 

An explicit solution for the diagonalization in case of 
k = 2 is well known in terms of the Clebsch-Gordan 
coefficients, which can be iteratively extended to higher 
k. We obtain 



Vj= Vj, 



£ s 7Ti£ rn s 



(12) 



£ s m£ m s 
\rrh£\<£,\m a \<s 



such that for all tp € Vi. 



£ s mp m s i 



L 2 ip = l(£+l)ip, 
S 2 ip = s(s + 1) tp, 



L z ip = nig ip 
S z ip = m s (p. 



Note that Vi,i S mim 3 may be zero for some £, s,me,m s 
Assembling equations ([TO and (fT2|) . we obtain 



£ni,...,n fc 

s rat m s 



£ s me m s 
• = CO Vl,£ s mi m s ■ 



That is, we have decomposed the configurations into the 
simultaneous eigenspaces of the angular momentum and 
spin operators ([5]). 



C. Restriction to fixed me and m s 



where TlV) <xl an< ^ -'-"iV'Xxl are the one " anc ^ two-body re- 
duced density matrices of the TV-body matrix (x\, re- 
spectively. Here ho is the single-particle (hydrogen-like) 
Hamiltonian and v ee is the interelectronic Coulomb po- 
tential, 



JAx-H' V ee = | 1 | ■ (14) 

2 \x\ \x-y\ 



Since these operators are independent of spin, we may 
effectively "trace out" the spin. With the standard no- 
tation 

(ab\cd) :— [ a(x\)b(xi) - c(x2)d(x2) dxiX2, 

Jrb \X 1 -X 2 \ 

(15) 

we obtain 

(X | H | ip) = tr h j m{xl +tr v ee T m{xl , (16) 



with 



(ho) ~{i\h \j), (17) 
(7|VXxl)i,- 1 7|i»<x| li a )> ( 18 ) 



(19) 

(t W{x \) := 51 (j«»^l r W(x|l <a >*^)- ( 2 °) 

ia<fc/3 



From general results about the angular momentum and 
spin algebra, it is well known that within an irreducible 
Z/ 2 -S l2 -eigenspace, the ladder operators L± = L x ± iL y 
and S± = S x ±iS y traverse the L z and S z eigenstates, re- 
spectively. Additionally, the ladder operators commute 
with the Hamiltonian H in ([3]) as well as with the CI 
Hamiltonian. Thus, in terms of eigenvalue determina- 
tion, it suffices to restrict to LS eigenstates with fixed 
m@ and m s . We adopt the convention in Ref. [iH . and set 
mg = 0, m s = s in the sequel. 

D. Reduced density matrices (RDMs) 

In this subsection, we will incorporate RDMs (see e.g. 
Ref. H and and |21| ) into the algorithmic framework to 
gain computational speedups and memory storage sav- 
ings. In fact, we use RDM's of wavefunction pairs. 

For any pair of states tp and X in the -ZV-body Hilbcrt 
space (Q]) , the matrix element of the Hamiltonian ([3]) can 
be rewritten as 

(x\Hip) = ti<u [h j W ( x \] +tr A 2 W [v ee r^)< x |] , (13) 



Here i, j, k, £ denote spatial orbitals and a, /?, 7, S are as- 
sociated spin-parts. The inequality constraint in the last 
sum refers to lexicographical ordering of spin-orbitals. 

By choosing the spatial orbitals real-valued, it fol- 
lows that (ij I hi) = (ji | kl) and (ij | ki) = (ij | Ik) for 
all i,j,k,£. Thus, together with (ij \ k£) = (k£\ij), 
it suffices to calculate (ij \ k£) for i < j, k < £ and 
(i,j) < (fcj^) (in lexicographical order) only. 

For our purposes, the following two features of the 
above RDM formalism are crucial. First, it avoids having 
to set up the full TV-particle operators Ho and V ee , allow- 
ing one to work instead with the one- and two-particle 
operators ho and v ee ; this leads to significant storage sav- 
ings, see Section [VI Dl Second, the map from ip and x to 
ri^w x | is an algebraic coefficient mapping which only de- 
pends on the symmetry types of the orbitals (i.e., s, p, d, 
...) and neither the radial wavefunctions nor the dilation 
parameters Z n g. So the F\-if>)( x \ can be precomputed for 
each angular momentum and spin symmetry eigenspace, 
without any reference to the Hamiltonian. The dilation 
parameters only enter the stage via the Couloumb inte- 
grals in v ee . 
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IV. HANDLING SLATER ORBITALS 
(STO'S) 

A. Ort honor malizat ion 

In this subsection we formalize the orthonormalization 
calculations for Slater-type orbitals (STOs) employed in 
Ref. [l6l equation (31). There, only Is, 2s and 2p wave- 
functions are considered, whereas here, we handle arbi- 
trary subshells. 

More concretely, the wave functions are given by 

/n-i-i \ ^ 

1pn£m(x) = SrU T l Y lm (d, If) ^ Kl,i C n £,i T l e %*~ T , 

r = \x\, £ = 0, 



,n — l, n = 1, 2, . 



(21) 



with b n £j being the ith coefficient of the associated La- 
guerre polynomial p(r) = I?^l_ x (^f), 

n + £ \(-2/n) 1 
2£+l + i) i\ 

and to-be determined orthogonalization coefficients 
c n i,i G K (i = 0, . . . , n — £ — 1) as well as orthonormaliza- 
tion constants s n i > 0. Since the spherical harmomics 
Ye m are orthogonal, we may fix the angular momen- 



tum quantum numbers £, m. Now using J c 
A "} t , orthogonality translates to 

= (ipnim I Ipkim) 



dr 



(x)ipk£m(x)r sinz9d??d<pdr 
(i+j + 2£ + 2)l 



s n esk£ } , o ne ^o k ijc n e,iC k ej . i+i+2i+3 



%iSm {c n i | B n iH^ k Bke c k t) 



(22) 



for all k = £ + 1, . . . , n — 1 . Here we have extended the 
vectors (cu,i)i=o,...,k-i-x by c k t,i = for i > k - £. The 
Hankel matrix H^ k is defined by 



H nk '■= (af +J (A)) . 



a\(X) := 



(2 + 2(1+1))! 

^4+2^+3 



and B n £ is the diagonal matrix diag (b n g^) i . Summarizing 
Eq. ([22]). we obtain 



cw JL Span (B nl B l nk B kt c ke ). 



k=l+l,...,n-l 



(23) 



so c„£ can be calculated iteratively for n = £+1, £+2, . . . , 
starting from the convention c^ + iy _ = 1. 

Note that the a\ are the moments of a nonnegative 
measure m on the positive real axis K+. Namely, let 
dm x ,i(t) = t 2 ^ +1 1 e" At dt, then 



t dm x ,i(t). 



The Stieltjes moment problem^ states that this is equiv- 
alent to the quadratic form given by H^ lk being positive. 

Once all c n i have been obtained, we may plug k = n 
into (|22p to calculate the normalization factors s n i from 



1 = (ipnim I Ipnim) = \ s nl\ 2 (c n £ | B n iH^ n B nl C n l) ■ (24) 



B. One-body integrals 



The one-body matrix elements 

(ipnim I tpn'i'm') '■= (ipnim I h \ Ipn'l'm') 



TrVt/Wm ■ ^ Ipn'l'm' ~ J—\ VWm tpn'i'm' d X 
2 \x\ ) 

(25) 



can be evaluated symbolically from Def. ([21]) by a com- 
puter algebra system, via symbolic differentiation and 
exact integration in spherical polar coordinates. 



C. Two-body integrals 

1. Switching to real-valued, cartesian 
coordinates 



As mentioned in Section [ill Dl we save computational 
costs by switching to real-valued spatial orbitals when 
calculating Coulomb integrals. Thus, for each fixed £, 
we apply a unitary base change Ui — (ui^ mm i) mm > to the 
spherical harmonics Ye m of degree £ to obtain real- valued 
polynomials in Cartesian coordinates, 

Zt m (x) := Ut, mm 'Ylm' 
m'=i,i-l,...,-i 

where C£ m>p G R, x q := Yl i=i xf. Pluggin this into 
Eq. ([2T[) results in real-valued Slater-type orbitals given 
by 



(n-l-l \ 



(26) 



where we have set d n n^ := bni,i Cnt,% to shorten notation. 

Concretely, for I = 1 we adapt Ref. UM and choose (in 
this order) 

1 [3 

(Z lm (x)) m = (pz,px,py) := -J- (x 3 ,x 1 ,x 2 ) . 



For £ = 2, 



(Z 2m (x)) m = (dO, dz, dm, dx, dy) : 



1 /15 



4 V 7T 

, 2xix 2l x\ — x 2 , 2x 2 x 3 , 2xix 3 



zx 3 x l J, 2 o 2 2 
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or — figuratively dx ~ 2 y z, dy ~ 2xz, dz ~ 2 x y, dO ~ 
(3z 2 — r 2 ) /v3j and dm ~ x 2 ~ y 2 . The corresponding 
unitary Ug read 

/ \ 1 / o o o o' 

rr 1 / o V20 \ 1 I -i i 

Ui = — = -101 and (72 = — f= 10001 
V2 V i i/ V^loioio 

\ -1 1 0, 

when arranging the spherical harmonics Yi m with de- 
creasing quantum number m. 



we obtain for the spinless Coulomb integrals (fT9|) with 
orbitals §n§ 



(Vee) ij>M = (Wfij I <Pk<Pt) = EE Cv ' 



(29) 



with Cy tqi A as in Eq. (|27[) and <v,g, A' the analogous 
constants for (pi~(x)ipt(x). 



2. Transformation to Fourier space 

We adapt the idea in Ref. [l6| to calculate Coulomb 
integrals between pairs of spatial orbitals by applying 
Fourier transformation. We use the normalization-factor- 
free convention 

GF/)(fc):= f f(x)e- ik -dx. 

Given one-electron orbitals y>i,<P2, ■■■ with tp i and Tifi € 
L 2 (R 3 ) n L°°(R 3 ), let f(x) := <pi(x)<pj(x) and g(x) := 
ifk(x) <pg(x). Then (see e.g. Ref. 1 161) 

( TO 1 <p m ) = -L / -L(j-/)( fe ) (fe) d 3 fc. 



3. Application to dilated Slater-type orbitals 

Taking pairwise products of the wavefunctions (1261) in- 
volves the convolution of coefficients, 



f(x) := ^nimCaO^n'^m'O 15 ) 

= s„£S„'£' Zi m {x)Zt> m i{x) 



r e 



with the discrete convolution 



i—k- 



Since we have switched to real- valued Car tesian orbitals Similar reasoning applies to the product Zt m Z e > m >, 



in the previous subsection, f(x) = tpi (x) ifj (x) can be 
expanded as 



Zt m (x)Z i/m ,{x) = ^2 ( c «m * Q'm')p ' xV ■ 
\p\ 1= t+t' 



!/=0 \9l ,92,93=0 



c„. z q e- Ar , r = \x\ 



Let 



(27) 

with constants c v ^ q and A > 0. Directly from the defini- 
tion of the Fourier transformation, it follows that 



g(x) := ^Mrh( x )' t Pn'i'm>( X ) 



(J-/) (fc) = (-1) 



be another pairwise product of wavefunctions. Then the 
Coulomb integral of these pairs equals 



" i qi+q2+q3 ^tt7 (J"e" Ar ) (fe), 



dk q 



(28) 



x - y\ 



where we have used the notation 



d q d qi 

7^ := n^T for each q G Nq. 



dk 

i—l 

It is well known that 

-Xr\ 



E ( 

\p\i=*+t' 



)p E 



(J"e~ Ar ) (fe) 



8Att 



(A 2 + P)' 



r, fc=|fe| 



(9 l 83 



dA* 



(30) 



Thus, precomputing the following integral over polar co- 
ordinates 



I q , q >(\,\') := (-i)<? 1+ « 2 +<» i«i+«i+«£- o 2 



1 

2~7 2 



V. COMPUTING THE CI LEVELS AND 
STATES 



00 ,tt , 27r / Qq 

n Jo 7o V^fc 9 (A 2 +P 
<9 q ' 8AV 



dfc 9 (A' 2 + fc 2 )' 



sin $ dip d$ dk 



Our overall algorithm for the CI method in Section lll CI 
consists of a symbolic part, symmetry reduction and re- 
duction to two-body space, and a numerical part, Hamil- 
tonian matrix diagonalization and orbital exponent opti- 
mization. 
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A. Symbolic precomputation 

The following pre-computational steps will allow us to 
calculate the matrix representation of the Hamiltonian 
quickly, given plug-in values for the dilation parameters 

1. Compute the simultaneous eigenspaces of the op- 
erators (0, via the algorithm described in Sec- 
tion nm 

2. For any simultaneous eigenspace of the opera- 
tors ([5|), choose an orthonormal basis (ipi, ...,ip r ) 
and calculate the one- and two-particle reduced 
density matrices TlViX^I and r IV>.><0 3 h respec- 
tively, of the iV-particle states \tpi) ('ipj\ for all 
i,j — l,...,r. Subsequently, trace out the spin 
part as defined in (|2"0"f and ([T5]l to obtain 7|^ 4 )(^.| 

and r | V» > | - 

3. For the Slater orbitals (l2"TT) . calculate symbolic ver- 
sions of the orthonormalization constants c n e and 
s n g in Section IIV Al via equations (|2"3")l and (|2~41 , 
respectively. Note that these constants still de- 
pend on the dilation parameters Z n i, which will 
be plugged in at the numerical optimization step 
below. 

4. Calculate symbolic matrix representations (fl~7f and 
(IT51) of the single-particle and electron-interaction 
Hamiltonians ho and v ee , using a computer algebra 
system and Eq. (|23|) for ho and Eq. for v ee . 
These matrices still depend on the orthonormaliza- 
tion constants s n g and c n g^ from Step 2, and on the 
dilation parameters Z n g. 

B. Numerical diagonalization and energy 
minimization 

For any given set of orbital exponents Z n e, we can now 
calculate and diagonalize the matrix representation of the 
Hamiltonian projected onto any LS-eigenspace, by using 
the reduced density matrix formalism in section IIIID1 In 
mathematical terms, 

1. For a current numerical value of the orbital expo- 
nents Zi t o, Z2,o, Za,i, evaluate the symbolic or- 
thonormalization constants s n g and c n g^ and the 

symbolic matrix elements of ho and v ee . 

2. Equation (flH)|) yields the matrix elements of the 
Hamiltonian on an LS-eigenspace with orthonor- 
mal basis (■01 , . . . , if) r ) , namely, 

(ipi\H\ipj) =tr (M^Wd) + tr (' 3 eef , |V J ><Vil) • 

(Note that it would be theoretically possible but 
computationally inefficient to carry out this step 
symbolically.) 



3. Obtain the ground state energy E = 
A min ((^|# l^) lJ= i,..., r ) 

4. Iteratively repeat these steps for different values 
of the orbital exponents within a suitable opti- 
mization routine to minimize the ground state en- 
ergy numerically. (We used a gradient-free simplex 
search method.) 



VI. COST ANALYSIS 

In what follows we review the computational speedup 
of the central algorithmic steps as compared to operating 
directly on the full TV-particle Hilbert space f\ N T-L. 



A. Configurations 

In this paragraph, we quantify the savings by the 
configuration calculus introduced in Section IIII Al To 
shorten notation, set gj := dimy„ ., and assume that the 
total particle number N is fixed. Thus, the dimension of 
the full iV-particle Hilbert space equals (^f 3 )- Consider 
configurations C ni,, " ,Tli with J2 n j = N. They partition 
the Hilbert space, and accordingly 

e dim ( c--)= £ n© 

m,...,n k m,...,rik j V 3/ 

= ( E /)=dim(A^) 

as expected. Now, assume we are given an algorithm of 
order C(dim p ), like, e.g., LS diagonalization with p = 3. 
Running this algorithm either applied to all configura- 
tions separately or to the full iV-particle Hilbert space 
incurs computational costs of order 

^2 dim(C ni '- ,n *f as compared to dim(A N H) P 

m,—,n k 
J2n j= N 

. ( 31 ) 

In what follows, we derive an estimate of the quotient of 
these two terms. The Stirling approximation of factorials 
and a logarithmic series expansion leads to 

f o\ ^ 1 j-r, , n 1 (g- 2 ") 2 
i y \ K2-- +9 (irg)-2e 

Plugging this into the left hand side of (f3~Tj) yields 

e n(*H"7°>-5>) 

m,...,n h j \'h/ J J-oc x ' 

x JJ^Ci+w) (7r 5j ) _f e~ pl "'* a S dm--- dn k . 
3 
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The Fourier transform of these integrals is the pointwise 
product of the individual Fourier transforms. One ob- 
tains 



2s e 



-int 



dn 



= (2 /tt) 5 (p_1) p- 5 2 P5 g3 e -9 i ( 4i P+*)/( 8 P) 

for each individual transform. Now, the inverse Fourier 
transform of the pointwise products gives the desired ap- 
proximation of the left hand side pip, namely 



(2/7r) i(i-*+M p §(i-*) (5prod) 



l(i-p) 



x 2 P9 — ( 5sum )- 5 e- 



-2N) Z 



where we have set g s 



T,j9j and 5pr 



(32) 



to shorten notation. Finally, dividing the (Stirling ap- 
proximated) right hand side of (j3"Tj) by (j3"2"]) yields the 
sought-after quotient 



(I 



i(fc-i)( P -i) 



-1) / ffprod \ 
\ Ssum / 



Note that this factor is independent of the particle num- 
ber iV. It equals 1 for p = 1, as expected. 

As concrete example, consider Chromium with three 
active subshells 3p, 3d, 4s, i.e., all subshells up to 3s are 
completely filled. Thus, in terms of the computation pa- 
rameters we have (51,32,53) = (dirnVp, dimVd, dimVg) = 
(6, 10,2), N e s = 12 (electron number in active orbitals) 
and algorithmic order p = 3, say. Then, the approxi- 
mated quotient equals 5ir 2 = 49.348, which is close to 
the exact number 50774322144/938076521 = 54.1. 



assume 0(1) cost for each of these Kronecker products. 
Summing up (13"3")) for all £ yields 



^ numcc {V12,, 
t=\lx-li\ 



(2h + 1) (2£ 2 + 1) 



(2*2 + 1) 



If 



(34) 



< (2£ 2 + 1) dim (V t ® V 2 ) . 



Now consider irreducible angular momentum eigenspaces 
Vj, j — 1, . . . , k, with respective quantum numbers £j. 
The computational cost of the iterated Clebsch-Gordan 
method will be dominated by the calculation of the total 
angular momentum eigenspaces of 



2li + 1, 



(01; , ::, ) ^ dimjV, , : , ) 

where each V 1 fe _ 1 .|. is an irreducible angular momem- 

tum eigenspace in Vj such that Ylift&i + 1) = 
dim(Vi<g)- • -<E> Vit_i). According to (f3"4)) . this requires not 
more than (2£fc + 1) • dim(Vi ® • • • ® Vfc) Clebsch-Gordan 
coefficients and associated Kronecker products. Thus, in 
case of all Vj being of uniformly bounded dimension, i.e., 



ii,---,*k < W, the cost is of order 
costcc (Vi, . . . , V k ) = (dim (Vl « 



^4)) • (35) 



So indeed, the cost is (almost) of the order of the problem 
size. 

The analysis for spin states is exactly the same, and 
the angular momentum and spin operators can be treated 
independently. Thus, the result (|35[) remains valid when 
considering both angular momentum and spin. 



B. LS diagonalization for sparse vectors 

In this paragraph, we show that the cost of the decom- 
position (fT2)) essentially scales linearly in the problem 
size dim(Vr), assuming a sparse structure of the associ- 
ated coefficient vectors. 

First, consider two irreducible angular momentum 
eigenspaces V\ and V 2 with quantum numbers £j and 
dimensions (21 j + 1) (j = 1, 2, without loss of generality 
h > £2)- Then the Clebsch-Gordan method partitions 
V\ ® V2 into total angular momentum eigenstates, i.e., 

V 1 ®V 2 = V X 2, h dim (Vi2,i) = 2£ + l. 
e=\ii-i 2 \ 

Each Vi2,e requires the computation of exactly 

numcG (Vi2,i) , 

(33) 

= (l x + i 2 + 1) (2£ + 1) - (£ x - i 2 f - £(l + 1) 

Clebsch-Gordan coefficients and Kronecker products (pi® 
f2 (fj ^ Vj). Due to the mentioned sparse structure, we 



C. Diagonalization of the Hamiltonian 

We now consider the exact reduction steps introduced 
in subsections III Al and IIII Cl the Hamiltonian can be di- 
agonalized within each LS eigenstate separately, and only 
states with quantum numbers mi = 0, m s = s need to 
be taken into account. (Partitioning into configurations 
is advantageous for the LS diagonalization only, since 
the Hamiltonian mixes configurations.) The latter saves 
a factor of (21 + 1) • (2s + 1) states with each L 2 -S 2 
eigenspace. The former, in the examples in Section IVIIl 
reduces the number of states by a factor of 10 2 to 10 3 . 

We illustrate the huge cost reduction by the ex- 
ample of the Chromium 7 S states with configurations 
[Ar] 3dMs 4p k 4d e such that j + k + I = 5 (see Sec- 
tion lVIip . The dimension of the full CI state space equals 
(5) = 65780 (since 5 valence electrons have to be al- 
located to 10+6+10 possible orbitals). By contrast, re- 
stricting to a typical symmetry subspace of interest, such 
as 1 S (i.e., L 2 — and S 2 — 3(3 + 1)), reduces the dimen- 
sion to 98, and taking S z maximal and L z — reduces it 
further to 14. 



11 



D. Storing RDMs instead of N-particle 
wavefunctions 

Even though the number of required wavefunctions has 
been reduced, each individual iV-electron wavefunction 
on a iCorbital space still requires, a priori, (^) entries. 

First — as illustrated in Section MI Dl — this cost 
can be reduced since the components of the Hamilto- 
nian matrix on a given iV-particle subspace only re- 
quires knowledge of the two-particle density matrices of 
any pair of iV-particle basis functions. These RDMs 

have (^) = 0(K A ) entries, namely (^\ij>}(x\)ij,kl with 
1 < i < j < K and 1 < k < t < K. 

Second, applying the spinless density matrix defined 
in equation (|2"0|) reduces the number K of single-particle 
orbitals by one half. 

Third, we note that the density matrix typically ex- 
hibits a sparse structure, so we actually need far fewer 
entries. This is related to prior LS diagonalization on 
the iV-particle Hilbert space. More precisely, the two- 
particle RDM of an iV-particle L 2 -S 2 -L z -S z eigenstate 
must commute with these symmetry operators on the 
two-particle space. Reconsider, for instance, the S 
states of Chromium with configuration [Ar] 3d 3 As Ap k Ad 1 , 
j + k + I = 5. A general spinless RDM with orbitals 

up to 4d has ( 23 ^ hl ) = 76176 entries. By contrast, the 
14 2 = 196 RDM's of the 14 7 S states with S z maxi- 
mal and L z = turn out to have, on average, only 94.3 
nonzero entries, the maximum number of nonzero entries 
which occurs being 648. 



VII. ANOMALOUS FILLING OF 4S AND 3D 
ORBITALS IN TRANSITION METAL ATOMS 

We have applied the algorithmic framework reported 
above to the calculation of ground and excited states and 
levels in 3d transition metal atoms. These continue to 
offer substantial computational challenges, due to the ir- 
regular filling of 4s versus 3d orbitals, strong correlations, 
and non-negligible relativistic effects. 

Previous computations have led to different results, 
depending on the level of theory used. Limitations of 
single-determinant Hartree-Fock theory for these atoms 
are discussed in Ref. [U Multi-determinant HartreeFock 
(HF) energies for the experimental ground state con- 
figurations (but not for competing configurations) are 
given in Ref. (with the exception of Cr), Ref. [24] 
(onl y fo r atoms with anomalous filling such as Cr), and 
Ref. [23. The interconfigurational ordering of 4s 1 3d" ver- 
sus 4s 2 3d n ~ 1 is discussed in Ref. UH for relativistic HF 
and in Ref. [ID and [l4| for DFT. Among the transition 
metal series Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, rela- 
tivistic HF rendered 4s 1 stable for Cr, Mn, Fe, Ni, Cu, 
even though experimentally only Cr and Cu have a 4s 1 



Atom 


Sym 


dim 




energy [a.u. 


1 

\ 




CI 


exp 


(subsp) 


01 


exp 


IVlUrlr 


K 


2 S 


2 S 


1 


-596.7993 


-601.9337 


-599.16478 


Ca 


x s 




2 


-674.2442 


-680.1920 


-676.75818 


Sc 


2 D 


2 D 


4 


-756.8908 


-763.8673 


-759.73571 


Ti 


3 F 




5 


-845.1599 


-853.3503 


-848.40599 


V 


4 F 


4 F 


4 


-939.1657 


-948.8394 


-942.88433 


Cr 


5 D 


7 S 


3 


-1039.0409 


-1050.4914 


-1043.3563 


Mn 


6 S 


6 5 


1 


-1144.9715 


-1158.2670 


-1149.8662 


Fe 


5 D 


3 D 


1 


-1256.7813 


-1271.6930 


-1262.4436 


Co 


4 F 


4 F 


2 


-1374.8903 


-1393.3526 


-1381.4145 


Ni 


3 F 


3 F 


1 


-1499.3759 


-1520.6907 


-1506.8709 


Cu 


2 D 


2 S 


1 


-1630.3692 


-1655.1317 


-1638.9637 


Zn 


4 S 


4 S 


1 


-1768.0729 




-1777.8481 



TABLE II. Ground state symmetries and energies from K 
to Zn predicted by minimal asymptotics-based CI (this pa- 
per) with active space [Ar] 3<f 4s fe and compared to experi- 
mental data^. Boldface denotes deviation from experiment. 
The dimension of the joint eigenspace of the symmetry op- 
erators (JS| which contains the unique ground state with 
L z = 0, S z maximal is denoted by dim. Also shown are 
multi-determinant Hartree-Fock energies for the experimen- 
tal grou nd state symmetries^ (for Ti and Cr see also Ref. l24l 
andlH). 

ground stateQ DFT does not fare better, regardless of 
the type of exchange-correlation functional used: 4s 1 is 
rendered stable by Becke 88 for Ti, V, Cr, Ni, CuM, by 
the local density approximation and Perdew-Wang for V, 
Cr, Co, Ni, Cu 13 ' 14 , and by B3LYP for V, Cr, Co, Ni, 
Cu 1 ^. The poor atomization energies of DFT functional 
such as Becke 88 and B3LYP for transition metal dimers 
(those for &2 even come out with the wrong sign) have 
been associated-^ to poor interconfigurational energies of 
the atoms. It is then of interest to revisit the latter from 
alternative theoretical points of view. 

Our results for the asymptotics-based CI model (A), 
(B), (C) in Section Hi CI are as follows. First, we consid- 
ered a minimal model for the third period elements K to 
Zn with configurations [Ar]3d 3 4s k , that is to say in the 
language of Section III CI we choose the cutoffs 

(n, *)mta = (3, 1) = 3p, (n, ^) max = (4, 0) - 4s. 



In fact, for Ni the experimental classification as 4s 2 should be 
viewed with some caution. A look at the actual dataSS shows 
that for Ni, relativistic J splittings are of the same order as 
the interconfigurational gap, and while the experimental ground 
state is a particular J state of the 4s 2 ( 3 -F) configuration, 4s 1 
( 3 D) would become stable if one averages over J according to 
multiplicity. 
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It turns out that the ground states from Ca to Zn always 
put two electrons in the 4s subshell, i.e. have configura- 
ton [.Ar]3eP4s 2 . See Table Ull Thus minimal asymptotics- 
based CI coincides with the empirical Madelung rule 
(which states that the subshells are filled in the order 
of increasing n + £ and, for equal n + £, in the order of de- 
creasing n) . Experimentally, this means that the method 
fails for the two anomalous atoms Cr and Cu. 

Next, to address this issue we enlarged the CI subspace 
for the series Ca, Sc, Ti, V, Cr by the higher subshells 
Ap and Ad. That is to say we changed the cutoffs to 

(n, £) min = (3, 1) = 3p, (n, £) max = (4, 2) = id, 

hence including all configurations [Ar]4d J As k Ap 1 Ad m , and 
restricted to k = 1 (4s 1 ) and k = 2 (4s 2 ), respectively. 
In each case, we considered only the L and S values se- 
lected by Hund's rules (i.e., we minimized first S and 
then L, taking into account one s and d subshell as in the 
minimal model above) , computed the corresponding sym- 
metry subspaces via the algorithm in Section fill Bl and 
determined the associated eigenstates and energy levels. 
The results are shown in Table IIIII 

Despite the smallness of the radial basis set, the pre- 
dicted ground state configurations and spin and angular 
momentum quantum numbers are in full agreement with 
the experimental data. Physically, interesting insights 
can be gained from the orbital exponents in Table IIII) 
and from the coefficients of the different configurations 
contained in the ground state. First, for Ca, the 4s elec- 
tron is more tightly bound than any d electrons, whereas 
for Sc, Ti, V, Cr, this effect is reversed, in both the 4s 1 
and the 4s 2 configuration, with 4s outside of both 3d 
and Ad. Second, considering for instance the 4s 1 ( 7 S) Cr 
ground state, the configurations and weight coefficients 
of the fourteen contributing basis states spanning the S, 



rriL = 0, ms = 3 


symmetry subspace of 3d 3 4s 1 Ap Ad' 1 


are as follows: 




3d 5 4p°4d° 


0.36 


3d 4 4p°4d 1 


0.63 


3d 3 4p 2 4d° 


0.056 


3d 3 4p°4d 2 (2D) 


0.31 and 0.50 


3d 2 4p 2 4d 1 (2D) 


0.036 and 0.038 


3d 2 4p°4d 3 (2D) 


0.17 and 0.28 


3d 1 4p 2 4d 2 (2D) 


0.016 and 0.014 


Zd}Ap°Ad A 


0.096 


3d°4p 2 4d 3 


0.0036 


3d°4p°4d 5 


0.012 



In particular, no configuration dominates, and the high- 
est weight configuration is not the naively expected 3d 5 
which one would enforce in both single-determinant HF 
and (L-S-adapted) multi-determinant HF, but 3d* Ad 1 
(weight 0.63), followed by 3d 3 Ad 2 (0.59), 3d 5 (0.36), and 
3d 2 Ad 3 (0.33). The highest- weight Cr 7 S basis function 
in which one of the 3d electrons has migrated to a Ad 



orbital is 

( \3d23di3do3dniAsAdra) — \3d 2 3di3do3d r ^AsAd rll ) 
+ \3d23d\3dni3draAsAdQ) — |3d23<io3d„i3<i I ,24s4(ii) 

+ |3di3d 3d nL 3d rf2 4s4d 2 ))/V5, 

with expressions of similar type for the remaining 13 basis 
functions. Despite the simple treatment of radial orbitals 
here, our results provide clear evidence of strong 3d-Ad 
inter-shell correlations in Cr, and suggests (by comparing 
energies of Tables ITTlandl lII|) a huge, symmetry-reversing, 
correlation energy in Cr of the order of 1 a.u. 

For more quantitative conclusions the radial basis set 
used here is too small, as is illustrated by our sys- 
tematically higher energies compared to the large-basis 
MDHF energies in Table [TTJ Our results constitute, how- 
ever, an important step towards an accurate quantita- 
tive computation of the correlation energy. The remain- 
ing step, which lies beyond the scope of the present 
paper, is to combine the exact lowest symmetry sub- 
spaces delivered by our symmetry reduction algorithm 
with high-accuracy, multi-parameter, self-consistent ra- 
dial orbital optimization routines as have been developed 
for Hartree-Fock theory^i 2 ^— . 

VIII. CONCLUSIONS 

We have developed and implemented an algorithm for 
CI calculations for atoms which allows full resolution of 
valence electron correlations in a large active space, via 
efficiently automated (and exact) symmetry reduction. 
Application to 3c? transition metal atoms shows that even 
very small radial basis sets yield the correct qualitative 
picture of the electronic structure when all correlations 
within and between the 3d, 4s, Ap, Ad shells are fully 
resolved and when orbital exponents are optimized self- 
consistently for the actual CI wavefunctions. We trace 
the qualitative accuracy of our results partly to the the- 
oretical fact that the asymptotics-based CI method used 
here yields the correct leading-order asymptotics for the 
low-lying spectral gaps in the fixed- N, large- Z limit. 

In subsequent work, we aim to obtain an accurate 
quantitative picture, by combining the careful algorith- 
mic treatment of correlations introduced here with suit- 
able large-parameter orbital optimization routines as are 
used in (numerical or Roothaan-type) Hartree-Fock the- 
ory 
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